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Extended free energy Lagrangians are proposed for first principles molecular dynamics simu- 
lations at finite electronic temperatures for plane-wave pseudopotential and local orbital density 
matrix based calculations. Thanks to the extended Lagrangian description the electronic degrees 
of freedom can be integrated by stable geometric schemes that conserve the free energy. For the 
local orbital representations both the nuclear and electronic forces have simple and numerically 
efficient expressions that are well suited for reduced complexity calculations. A rapidly converging 
recursive Fermi operator expansion method that does not require the calculation of eigenvalues and 
eigenfunctions for the construction of the fractionally occupied density matrix is discussed. An effi- 
cient expression for the Pulay force that is valid also for density matrices with fractional occupation 
occurring at finite electronic temperatures is also demonstrated. 

PACS numbers: 



I. INTRODUCTION 

Molecular dynamics simulations are widely used in ma- 
terials science, chemistry and molecular biology. Unfor- 
tunately, the most accurate molecular dynamics meth- 
ods that are based on self-consistent field (SCF) theory 
[T] are often limited by some fundamental shortcomings 
such as a very high computational cost or unbalanced 
phase space trajectories with unphysical hysteresis ef- 
fects, numerical instabilities and a systematic long-term 
energy drift [H Recently an extended Lagrangian 
framework for first principles molecular dynamics sim- 
ulations was proposed that overcomes most of the previ- 
ous problems [4] . The new framework was originally de- 
signed for ground state Born-Oppenheimer molecular dy- 
namics simulations based on Hartree-Fock [5^ or density 
functional theory [6] at zero electronic temperatures. 
In this paper we supplement previous work by propos- 
ing extended free energy Lagrangians for first principles 
molecular dynamics simulations also at finite electronic 
temperatures [8j|9|. 

Our extended Lagrangian free energy dynamics is for- 
mulated and demonstrated both for plane-wave pseu- 
dopotential and local orbital density matrix based cal- 
culations. For the density matrix formulation our ap- 
proach is given by a generalization of the original ex- 
tended Lagrangian framework [4] using a density ma- 
trix representation of the electronic degrees of freedom 
with fractional occupation of the states. The plane-wave 
formulation is based on the wave function extended La- 
grangian molecular dynamics method that recently was 
proposed by Steneteg et al. ^TU\ . 

First we formulate the extended Lagrangian approach 
to first principles molecular dynamics at finite electronic 



temperatures for local orbital density matrix representa- 
tions and thereafter for plane wave pseudopotential cal- 
culations. A numerically efficient expression for a gen- 
eralized Pulay force that does not require idempotent 
density matrices is demonstrated and a recursive Fermi 
operator expansion algorithm for the density matrix is 
discussed. Thereafter a Vcrlet integration scheme for the 
equations of motion is presented. Finally, we illustrate 
the extended Lagrangian free energy dynamics by some 
examples before giving a brief summary. 

II. FIRST PRINCIPLES MOLECULAR 
DYNAMICS AT FINITE ELECTRONIC 
TEMPERATURES 

The first principles molecular dynamics scheme pre- 
sented in this paper is a finite electronic temperature 
generalization of the extended Lagrangian framework of 
time-reversible Born-Oppenheimer molecular dynamics 
[31 H]. In many ways the generalization is straightfor- 
ward, except for an additional entropy term, problems 
arising in the construction of the density matrix, the cal- 
culation of the Pulay force for fractional occupations and 
the phase alignment problem in the integration of the 
wave functions |10j . Our two different formulations of 
extended Lagrangian free energy molecular dynamics are 
based on either an underlying local orbital representation 
or a plane wave pseudo potential description. Alterna- 
tive free energy formulations suitable for extended La- 
grangian self-consistent tight-binding molecular dynam- 
ics simulations [TTl [T^] will be presented elsewhere. 

A. Density matrix extended free energy 
Lagrangian 
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Within an atom-centered local orbital representation 
we describe our first principles molecular dynamics at 



2 



an electronic temperature, Te, by a density matrix (DM) 
extended free energy (XFE) Lagrangian, which we define 
as 



£gFf (R, R, P, F) = i Efc MkRl - C/[R; D] 
+ iMTr[p2] _ i^,,^2xr[(i^ _ p)2] + Te5[i?]. 



(1) 



The extended dynamical variables P and the velocity P 
in the Lagrangian are auxiliary electronic degrees 

of freedom evolving in a harmonic potential centered 
around the self-consistent free energy ground-state solu- 
tion D [4]. Here P, P, and the self-consistent free energy 
ground state D, are orthogonal density matrix represen- 
tations of the electronic degrees of freedom. The relation 
between D and the non-orthogonal (atomic orbital) rep- 
resentation, _D, is given by the congruence transformation 



D = ZDZ' 



(2) 



where Z is a congruence factor determined by the overlap 
matrix S (e.g. see Ref. [13 ) from the condition that 



SZ 



(3) 



The constants fi and lu are fictitious electron mass and 
frequency parameters. The potential U[R;D] is defined 
at the self-consistent field (SCF) electron ground state 
D, at which the free energy functional 



n[R,D] = [/[R; D] -T^S[Dl 



(4) 



has its minimum for a given nuclear configuration, R = 
{Ri}, [SI in] under the constraints of correct electron oc- 
cupation, Tr[D] = Nocc- The electronic temperature 
can be different from the ionic temperature Tion- We as- 
sume that U[Il; D] is the total electronic energy, includ- 
ing nuclear ion-ion repulsions, in self-consistent density 
functional theory, Hartree-Fock theory or some of their 
extensions that are based on an underlying SCF descrip- 
tion. The mean field entropy term in C^^^, 

S[D] ^ ~2kBTT[Dln{D) + (I - D)\ii{I - D)], (5) 

is included to provide the variationally correct energet- 
ics and dynamics [TH [15]. The factor 2 above is used 
for restricted calculations (not spin polarized) when each 
state is doubly occupied, which is assumed throughout 
the paper. 



1. Euler- Lagrange equations of motion 

The time evolution of the system described by C^^^ 
is determined by the Euler-Lagrange equations. 







Q/-XFE 


dt 


\ dRk J 


dRk 
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Qz-XFE 


dt 


K dP J 


dP 



(6) 
(7) 



which give the equations of motion for the nuclear and 
electronic degrees of freedom. 



dRk 



dRk 



dRk 



(8) 

fiP = ^Mj'^iD-P). (9) 
In the limit /i — > 0, we get the decoupled equations 



MkRk 



dU[Tl;D] , dS[D] 



dRk 



dRk 



P = uj^{D-P), 



(10) 
(11) 



where the self-consistent density matrix D is given from 
the Fermi operator expansion 



D = 



under the canonical condition of charge neutrality, 



(12) 



(13) 



i.e. the chemical potential fiQ is set such that the trace of 
D is the number of occupied states. The inverse temper- 
ature /3 = l/{kBTe). The orthogonafized efii'ective single- 
particle Hamiltonian H, i.e. the Fockian or Kohn-Sham 
Hamiltonian, is given by a congruence transformation of 
the Fockian in a non-orthogonal (atomic orbital) repre- 
sentation, H, where 



H = Z' HZ. 



(14) 



Since we use the limit /x — 0, we recover the regular 
free energy Lagrangain (without extention) and the total 
free energy is a constant of motion that is unaffected by 
the extended variables P and P. This is in contrast to 
extended Lagrangian Car-Parrinello type molecular dy- 
namics [TBHE]. The major advantage of the auxiliary set 
of variables P and P occurs in the SCF optimization and 
in the geometric integration of the equations of motion 
discussed in section III below. 



B. Wave Function Extended Free Energy 
Lagrangian 

In plane wave (PW) pseudopotential methods it is con- 
venient to propagate both the electronic wave functions, 
^ = {'/'m(r)}, and the density, p(r). In this case we de- 
fine the wave function (and density) extended free energy 
Lagrangian [TU] as 

/:^FE(R^ ^ <j,^ p^p) = \ J2k MkRl - C/[R; n-] 
Hf^Em I l0™(r)Pdr - J l^-(r) - 0„(r)pdr 



(15) 
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The auxihary electronic degrees of freedom $ and p are 
extended through harmonic oscillators centered around 
the self-consistent (sc) ground state wave functions ^^'^ — 
{^'m(r)} and densities rf^. As above, ^ is a fictitious 
electron mass parameter and a; is a frequency param- 
eter determining the curvature of the harmonic poten- 
tials. As for the density matrix extension, the potential 
[/[R; rf^] is defined at the self-consistent field (SCF) elec- 
tron ground state rf^ ^ at which the free energy functional 



17[R,n'''=] = C/[R;n'''] - Te5[f], 



(16) 



has its minimum for a given nuclear configuration, R = 
iSnSj under the constraints of correct electron occu- 
pation, / n*'^(r)dr = A"occ- The mean field entropy term 
in /3p™ is given analogous to Eq. ([s]) as a function of the 
occupation numbers, f = {/;}, where 

S[f\ = -2kB J2 (/^ Hfi) + (1 - /.) ln(l - /.)) , (17) 

i 

and 

= + (18) 

Here Ei are the eigenvalues of the effective single particle 
(Kohn-Sham) Hamiltonian, i.e. 



(19) 



The entropy term S[f] provides a variationally correct en- 
ergetics and dynamics [T31 [H] where the electron density 
is given by 

"''« = E//'"I^™WI'^^- (20) 



1. Euler-Lagrange equations of motion 

As for the density matrix extended free energy La- 
grangian, the dynamics for the wave function extended 
free energy Lagrangian is determined by the Euler- 
Lagrange equations, which give the equations of motion 
for the nuclear and electronic degrees of freedom, 

dU [R; n"'^] 




(21) 



2 dRk 

,2/ 



l\n''{r)-p{r)\^dr^ + 



dS[f] 



= - $), (22) 

Mr) = ^c.2(n-(r)-p(r)). (23) 

In the limit /i — > 0, we get the decoupled equations 



MkRk 



5C/[R;n-] dS[{] 



dRk 

* = w2(*"= - $), 
p(r) = c.2(n-(r)-p(r)). 



dRk 



(24) 

(25) 
(26) 



Also in this case, the total free energy is a constant of 
motion in the limit — 0. The advantages with the ex- 
tended equations of motion occur in the SCF optimiza- 
tion and geometric integration discussed in section III. 



C. Generalized Pulay force 



The nuclear force J"*"' = MkRk in Eqs. (10 1 and (24) 



can be partitioned into three separate terms. 



(27) 



The first term, , is the Hellmann-Feynman force term 
[20j l2Tj for a fractional occupation of the states [Ml [15] 
(including nuclear- nuclear repulsion). The second term, 
, corresponds to the Pulay force term 20J, and the 
third term, 



T^k^T,{dS/dRk), 



(28) 



is an additional entropy force term. 

In a plane wave formulation the Pulay force term usu- 
ally vanishes, but only at Tg = 0. At finite electronic 
temperatures the Pulay term is finite even for a constant 
plane wave basis set representation. In this case however, 
the Pulay term is exactly cancelled by the entropy force 
term J^^ if the entropy is given by Eq. (17). This is not 
true for a local orbital basis set representation. In this 
case only a part of the Pulay force term is canceled 
by the entropy contribution J^'^ f2?. Nevertheless, this 
cancellation provides a significant simplification of the 
Pulay term. The remaining force term can be viewed as 
a generalized (G) Pulay force [52], 



J-P + J-'S = 2Tr[S-'^HD{dS/dRk)], (29) 



which is valid for finite temperature ensembles with non- 
integer occupation [22J. For a constant orthonormal 
plane wave basis set representation the overlap matrix 
is equal to the identity matrix, i.e. S — I, and the 
generalized Pulay term disappears. A similar convenient 
expression for the Pulay force term was recently given 
by Schlegel et al. [TH], but without considering any en- 
tropy contribution from fractional occupation. A gen- 
eralized Pulay force was also discussed in connection to 
the efficient Ehrenfest-like molecular dynamics scheme 
by Jakowski and Morokuma [23j. In the present formu- 
lation the electronic entropy is necessary to make our 
free energy and forces variationally correct in a molecu- 
lar dynamics simulation with factional occupation of the 
states at finite electronic temperatures [TH [IS] . The den- 
sity matrix or the density that minimizes the free energy 
functional in Eq. (|4| or in Eq. ( 16 ) with the particular 
form of the entropy term in Eq. (|5|) or in Eq. (17) are 



given by the Fermi operator expansion in Eq. ( 12 ) or by 
Fermi weighted integration in Eq. (20). If other expres- 



sions of the entropy are used, the Pulay force expression 
in Eq. ( 29 ) may still be valid 24| , but the density matrix 
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D and wave functions will have a fractional occupation 
distribution of the eigenvalues that is different from the 
Fermi-Dirac function. 



III. GEOMETRIC INTEGRATION 

Thanks to the extended Lagrangian formulations of 
free energy molecular dynamics, Eqs. ([I]) and (15 1, it 



is possible to simultaneously integrate both the nuclear 
and the electronic degrees of freedom in Eqs. (10 1 and 



(111 for the density matrix formulation, or Eqs. ( 24 ) and 



(25) for the plane wave extension, with almost any kind 



of geometric integration scheme developed for celestial 
or classical molecular dynamics, e.g. the time-reversible 
Verlet algorithm or higher-order symplectic integration 
methods [H [TOj [25l [26^ Morevover, the auxiliary degrees 
of freedom P in Eq. (IT]) or $ and p in Eq. ( 15 ) will stay 
close to the free energy ground state D or ^''^^ and n^'', 
since P, $, and p evolve in harmonic potentials centered 
around D, 'i/^'^, and n^"^. It is therefore possible to use 
the extended dynamical variables as an efficient initial 
guess of D, or n^'^, in the iterative SCF optimization 
procedure that minimizes the free energy, i.e. 



D{t) ^ SCF [R;P(t)] 



{vI/-(i),n-(0} = SCF[R;<i>(i),p(i)] 



(31) 



This choice of initial guesses also avoids the fundamental 
problem of a broken time-reversal symmetry in the prop- 
agation of the underlying electronic degrees of freedom 
PHI], since the extended dynamical variables can be in- 
tegrated by a time-reversible algorithm, for example, the 
Verlet scheme [27]. This is in contrast to regular Born- 
Oppenheimer-like molecular dynamics, where the initial 
guess of the SCF optimization is based on some extrapo- 
lation of the self-consistent solutions from previous time 
steps [2, 28 -32 , for example 

D{t) = SCF [R; 2D{t - 6t) - D{t - 2St)] , (32) 

{■f"'{t),n"'{t)} = SCF [R; $(t - 6t), p{t - 6t)] , (33) 

which breaks time-reversal symmetry because of the in- 
complete and non-linear SCF optimization. A broken 
time-reversal symmetry leads to a hysteresis effect and a 
systematic drift in the total energy and the phase space 
[U 13]. Only by the costly procedure of increasing the 
SCF convergence to a very high degree of accuracy is it 
possible to reduce the drift, though it never really disap- 
pears. As will be shown below, our extended Lagrangian 
approach allows for highly efficient and stable free energy 
molecular dynamics simulations without any significant 
energy drift. Often only one single SCF iteration is suf- 
ficient to provide accurate trajectories. 

In principle, time-reversibility is not a necessary cri- 
terion for long-term energy conservation. A more fun- 



damental aspect is the conservation of geometric proper- 
ties of the exact fiow of the underlying dynamics as in 
symplectic integration schemes [H 133] • Here we only 
consider the Verlet integration algorithm, which in its ve- 
locity formulation is both time-reversible and symplectic. 
First we present the integration of the density matrix and 
thereafter the electron wave function and density integra- 
tion. 



A. Density matrix integration 

The time-reversible Verlet integration [57] of the den- 
sity matrix P{t) in Eq. (11), 



P{t + St) = 2P{t) - P{t ~ St) + 5t^P{t), 



(34) 



that includes a weak external dissipative force that re- 
moves the accumulation of numerical noise |34| , has the 
following form: 



Pn+l = 2P„ - P„_i + St'^LJ^iDn - Pn) 



K 

CfcPn-fc- 

fc=0 



(35) 

(30) In the first K + \ initial steps we set P„ = D„, where 
Dn = D{tQ + nSt), and we use a high degree of SCF 
convergence, typically 10-15 SCF cycles per time step. 



The last term on the right hand side of Eq. ( 35 ) is the 



dissipative force term, which is tuned by some small cou- 
pling constant a. This term corresponds to a weak cou- 
pling to an external system that removes numerical error 
fiuctuations, but without any significant modification of 
the microcanonical trajectories and the free energy [34j . 
A few optimized examples of the dimensionless parame- 
ter, St^uP , the Ck coefficients and a are given in Table I 
(see Ref. [31] for details). The integration coefficients 
in Table [I are optimized under the condition of stabil- 
ity under incomplete, approximate SCF convergence. A 
larger value of a in Table |l] gives more dissipation but also 
a larger perturbation of the molecular trajectories. For 
a = the scheme is exactly time-reversible and perfectly 
lossless without any dissipation of numerical errors. In an 
exactly time-reversible and lossless scheme small numer- 
ical errors may accumulate to large fluctuations during 
long simulations, which can lead to a loss of numerical 
stability. It is particularly important to include dissi- 
pation when the numerical noise is large, for example, in 
reduced complexity calculations utilizing an approximate 
sparse matrix algebra. 



B. Wave function and density integration 

The Verlet integration of the wave functions and the 
density in Eqs. (25) and (|26[), including a weak external 



dissipative force term that removes a systematic accu- 
mulation of numerical noise (as described above) has the 
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TABLE I: Coefficients for the Verlet integration sciieme with 
the external dissipative force term in Eq. (35 1. The values 
are derived in Ref. 34 , which contains a more complete set 
of coefficients. 



K 




ax 10"^ 


Co 


Cl 


C2 


C3 


C4 


C5 


C6 C7 





2.00 



















3 


1.69 


150 


-2 


3 





-1 








5 


1.82 


18 


-6 


14 


-8 


-3 


4 


-1 




7 


1.86 


1.6 


-36 


99 


-88 


11 


32 


-25 


8 -1 



following form: 



n+l 



Pn+1 



2$. 



K 



fe=0 



— mi 



K 

+ Oi'^CkPn-m- 
k=0 



(36) 



(37) 



Because of the cubic scaling of the computational cost as 
a function of system size, the diagonalization approach 
is ill-suited for reduced complexity calculations of large 
extended systems. In this case we instead propose us- 
ing recursive Fermi operator expansion methods, both at 
finite and zero electronic temperatures. In plane wave 
pseudopotential methods usually some form of iterative 
eigensolvers are used to calculate a subset of the lowest 
lying eigenf unctions |28) . Plane wave matrix represen- 
tations are not sparse and linear scaling of the cost as 
a function of the number of atoms can in general not 
be achieved. Because of the large number of plane wave 
basis functions, an explicit construction of the density 
matrix is difficult and it is only for local atom centered 
basis set representations for which it is useful to calcu- 
late the density matrix. Here we will discuss efficient 
recursive expansion methods for the construction of the 
density matrix at finite and zero electronic temperatures. 
The recursive Fermi operator expansion methods allows 
for very efficient reduced complexity calculations of large 
systems if sparse matrix algebra can be used. 



Sine the electronic ground state wave functions, are 
unique except with respect to their phase, a straightfor- 
ward Verlet integration of the wave functions may lead 
to large and uncontrolled errors. To avoid this problem 
we have included a subspace alignment through a unitary 
rotation transform [/„ in Eq. ( 36 1 |10j , where 



Un 



arg rmn||^'„C/' - $„||i., 



(38) 



where || • ||f denotes the Frobenius norm. The rotation 
Un that minimizes this distance between ^„f7' and $„ 
is given by 



[/ = (OOt)-i/2(9^ 



(39) 



where O = {^^^\^n) ^S\- Using this subspace alignment 
we can use the same coeflicients as in the density matrix 
integration that are given in Tab. [ij 



IV. RECURSIVE FERMI OPERATOR 
EXPANSION FOR THE DENSITY MATRIX 

The density matrix can be calculated from a diagonal- 
ization of the effective single-particle Hamiltonian, H, 
where 



(40) 



Here Q is an orthogonal matrix consisting of the Hamil- 
tonian eigenvectors and is a diagonal matrix with the 
corresponding Hamiltonian eigenvalues {si}. By replac- 
ing these eigenvalues {si} with the occupation numbers 
{fi} in Eq. (18), the density matrix is given by 



(41) 



A. Te > 

The density matrix at finite electronic temperatures 
used in the calculation of the nuclear and electronic 



forces, Eqs. ( 10 1 and (111, can be constructed by a recur- 
sive Fermi operator expansion [36], 



D = [e^i^-f^oD + jj ^ rmiv„-,-i{. ..r^"'\H) . . .)). 

(42) 

There are certainly numerous fairly well established al- 
ternative techniques that can incorporate the effects of 
the electronic distribution at finite temperatures, such 
as Matsubara, Green's function, integral representation, 
Chebyshev Fermi operator expansion, or continued frac- 
tion methods [57^4 7j . However, in contrast to schemes 
that are based on serial expansions, the recursive expan- 
sion in Eq. ( 42 ) leads to a very high-order approximation 



in only a few number of iteration steps. Here we will use 
a recursive Fermi operator expansion based on the Pade 
polynomial 



Vn{x) = X^lx^ + {1 - Xf] 



(43) 



in Eq. (|42]) for n > [SB]. 

The problem of finding the chemical potential fig in 
Eq. ( 42 ) such that the correct occupation is achieved. 



(44) 



can be solved by a Newton-Raphson optimization using 
the analytic derivative of the density matrix with respect 
to ^0, 



dD 

dpQ 



PD{I~D). 



(45) 
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The recursive Fermi operator expansion of the 
temperature-dependent density matrix, which automati- 
cally finds the chemical potential by a Newton-Raphson 
optimization, is given by Algorithm [T] The algorithm 
is a slight modification of the original scheme given in 
Rcf. [36 and was previously presented in Ref. [22]. A 
full derivation is given for the first time in Appendix [K\ 
The main difference to the original algorithm is that the 
expansion order m is independent of the temperature and 
that the expansion is performed for D directly, instead of 
the virtual projector I — D. Note also that the choice of 
m determines both the initialization of Xq = Vq"^^ (H) as 
well as the expansion order as is shown in the appendix. 
The initial guess for fiQ in the Newton-Raphson optimiza- 
tion has to be particularly good at low temperatures, i.e. 
when /? is large. If m is small the analytic derivative 
PD{I — D) is more approximate, which also may 



dP 
dfJ.0 



require a better initial guess. 



Algorithm 1 Recursive Fermi operator expansion. 



H < 

m < 

Mo 



- Z^HZ 

- Number of recursive iterations 
l/iksTe) 

- Initial guess of /io 
while Occupation Error > Tolerance do 

X^^\l-{H- ^o/)/?/22+™ 
for n = \ : m do 

solve [X2_i-h(J-X„_i)2]X„ = 
end for 

Occupation Error = |Tr(_D) — A'^occ 

MO ^ MO + [iVocc - Tr(Li)]/Tr[/3D(J - D)\ 

end while 

D ^ ZDZ'^ 



X„ 



|49) . shown in Algorithm [2j Here (t„ = ±1 and is chosen 
to minimize the occupation error, |Tr[X„+i] — Aocd- The 
idempotency error, 

Idempotency Error = |Tr[X„ -X^]] = |Tr[X„(/ - X„)] |, 

(46) 

which determines the degree of convergence, can be esti- 
mated from the size of the occupation correction, since 
|Tr[X„(/ - X„)]| - |Tr[X„+i] - Tr[X„]|, and towards 
convergence (e.g. n > 10) the continuous decrease of the 
idempotency error can be estimated from every second 
value. In the initialization, Emax and Emin are upper and 
lower estimates of the spectral bounds of H. To improve 
matrix sparsity, numerical thresholding can be applied 
after each matrix multiplication within controllable total 
error bounds ^13,, ,50r,52j . The recursive 2nd-order oc- 
cupation correcting expansion algorithm is fast, robust, 
and very memory efficient compared to other linear scal- 
ing methods ^3 " 56 . If more detailed information about 
the homo-lumo gap and the spectral bounds are avail- 
able, an efficient boosting technique based on shifting 
and rescaling of the second order projection polynomials 
in Eq. ( 46 1 can be applied to accelerate the idempotency 
convergence |57j . 



Algorithm 2 Recursive Fermi operator expansion at 
zero electronic temperature |49) 

H Z'^HZ 

Xo <— (Cmax — Emin) ^ (£maxl — H) 

while the Idempotency Error is decreasing do 

X„+i <— X„ + ari{X„ — X„) 
end while 

D ^ lim„_>oo ZX„Z^ 



In the recursive Fermi operator expansion. Algo- 
rithm [l] a system of linear equations is solved in each 
iteration. The numerical problem is well conditioned and 
the solution from the previous cycle Xn-i is typically 
close to the solution X„. A linear conjugate-gradient 
solver is therefore very efficient [3S1 [35] . The Fermi oper- 
ator expansion algorithm can be formulated based only 
on matrix-matrix operations and does not require a di- 
agonalization of the Fockian. It is therefore possible to 
reach a reduced complexity scaling of the computational 
cost as a function of system size if sparse matrix algebra 
can be utilized. 

B. Te = 

The recursive Fermi operator expansion. Algorithm [T] 
provides an efficient calculation of the density matrix at 
finite electronic temperatures, Te > 0. The algorithm 
can be applied also at zero temperatures, Tg = 0. How- 
ever, in this case there are more efficient recursive Fermi 
operator expansion schemes, for example, the second or- 
der occupation correcting spectral projection technique 



V. EXAMPLES 

A. Computational details 

The extended Lagrangian free energy molecular dy- 
namics schemes were implemented in FreeON, as suite of 
freely available programs for ab initio electronic structure 
calculations developed by Challacombe and co-workers 
[55] , and in VASP, the Vienna Ab initio Simulation Pack- 
age [5U] . 

For the calculation of the electronic entropy contribu- 
tion to the free energy we rely on a diagonalization of the 
Fockian or the Kohn-Sham Hamiltonian, which would 
not be possible within a reduced complexity scheme. A 
sufficiently accurate expansion of S in Eq. ([s]), which is 
suitable within a reduced complexity calculation may be 
possible, but it is currently an unsolved problem. For 
the construction of the density matrix we use either the 
recursive Fermi operator expansion. Algorithm [l] or, as 
a comparison, an "exact" calculation based on the con- 
ventional (non-recursive) exponential form of the Fermi- 
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FIG. 1: (Color online) The fluctuations in the total free 
energy, i?^* = i?K + U — TeS, of a water molecule using regu- 
lar free energy molecular dynamics for a conventional density 
matrix based Born-Oppenheimer method with a linear ex- 



trapolation of the electronic state, Eq. (321, in comparison to 



the density matrix extended Lagrangian free energy molecular 
dynamics. 



FIG. 2: (Color online) The fluctuations in the nuclear ki- 
netic and potential energy, Ek + U, in comparison to the 
total free energy, E^p^ = Ek + U — TeS, for a density matrix 
extended Lagrangian free energy molecular dynamics simula- 
tion of a single water molecule using Hartree-Fock (RHF) or 
hybrid density functional theory (PBEO i58l). The average 
nuclear temperature was approximately room temperature, 
i.e. Tion ~ 300 K. 



Dirac distribution in Eq. (42), which requires a full diag- 



onalization of the Hamiltonian. 

The plane wave pseudo potential calculations were car- 
ried out within the local density approximation (LDA) 
for the exchange correlation energy. Ultrasoft pseudo 
potentials were used [HI], and for the SCF optimization 
an iterative RMM-DIIS scheme [5H [53] was appHed. The 
plane wave cutoff was 246 eV and a mesh of 64 A;-points 
was used. The electronic temperature was set to = 500 
K around which also the nuclear temperature fluctuated. 



B. Free energy conservation 

1. Density matrix extension 

In a free energy simulation based on regular first prin- 
ciples molecular dynamics the self-consistent electronic 
state is propagated through an extrapolation from pre- 
vious time steps, which is used as an initial guess to the 
SCF optimization procedure as was shown in Eq. ( 32 1 . 



This approach leads to the breaking of time-reversal sym- 
metry with an unphysical systematic drift in the energy 
and phase space. Figure [T] illustrates this problem, which 
is avoided with the density matrix extended Lagrangian 
free energy molecular dynamics. 

Figure [2] shows the free energy conservation vs. the 
fluctuations of the sum of the nuclear kinetic and poten- 
tial energy, Ek + U. Only one or two SCF cycles per 
time step are used, which shows that the simulation is 
stable, conserving the free energy under incomplete ap- 
proximate SCF convergence. While the amplitude of the 



fluctuations in Ek + U is dependent on the electronic 
temperature Te, the amplitude of the total free energy 
oscillations, E*^*- ^ Ek + U - T^S, is not. This IS an 
agreement with the total free energy being a constant 
of motion as in the density matrix extended free energy 
Lagrangian, Eq. ([I]), in the limit /i — > 0. 

Figure |3] shows the corresponding behavior of the free 
energy fluctuations of the sum of the nuclear kinetic and 
potential energy, Ek + U for a single water molecule em- 
bedded in a CgQ Bucky ball. The electronic temperature 
was set to Te = 5, 000 K, while the nuclear temperature of 
the Bucky ball was initially zero, and the water molecule 
was given an initial velocity and hence had a non-zero ini- 
tial temperature. The time evolution of the total temper- 
ature, T, the temperature of only the Bucky ball, Tceo, 
and the temperature of only the water molecule, Twatcr, 
are shown in Figure |4j Clearly visible in the figure is the 
process of kinetic energy transfer from the water molecule 
to the Bucky ball; as the water molecule repeatedly col- 
lides with the wall of the Bucky ball, kinetic energy is 
transferred between the two sub-systems and the Bucky 
ball heats up, while the water molecule cools down. 

If we use the original expression of the Pulay term (50] , 



J-P = 2Tv[DFD{dSldRk% 



(47) 



which is valid only for an idempotent density matrix in 
the limit of Tg — 0, i.e. when D = DSD, we find that 
it affects the energetics even at fairly modest electronic 
temperatures. Figure [S] gives a comparison of the free en- 
ergy calculated either with the original zero-temperature 
expression in Eq. (|47[) or its finite temperature general- 
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FIG. 3: (Color online) The fluctuations in the nuclear kinetic 
and potential energy, Ek + U, in comparison to the total free 
energy, = Ek + U — T^S, for a density matrix extended 
Lagrangian free energy molecular dynamics simulation of a 
single water molecule embedded within a Ceo Bucky ball using 
Hartree-Fock theory. The average nuclear temperature was 
approximately Tion ~ 1000 K. 



FIG. 4: (Color online) The fluctuations in the total nuclear 
temperature, Ttotai, in comparison with the temperature of 
only the Bucky ball, Tceo and the temperature of the water 
molecule, T„ater- We start our simulation with the Bucky ball 
at Tceo = K and a water molecule at the center inside the 
Bucky ball with a non-zero initial velocity. While the total 
temperature remains fairly constant, the energy transfer from 
the water molecule to the Bucky ball is clearly visible in the 
rising Tbeo and the falling r„ater- 



ization in Eq. ( 29 ) . Only the generalized Pulay force has 



the correct variational properties derived from the ex- 
tended free energy Lagrangian, which gives a significant 
improvement in the conservation of the free energy. 



2. Plane wave extension 



Figure |6] demonstrates how the total free energy for 
the conventional Born-Oppenheimer molecular dynam- 
ics and the plane wave extended Lagrangian free en- 
ergy molecular dynamics behave over time in the plane 
wave formulation. A system of 8 silicon atoms are sim- 
ulated for a total of 5 ps using a time step of 0.25 
fs. The figure effectively depicts how the regular Born- 
Oppenheimer molecular dynamics exhibit an unphysical 
systematic drift in the total free energy. This issue is 
solved by the use of the extended Lagrangian free en- 
ergy formulation, which with equivalent settings other- 
wise shows no systematic drift. In Fig. [7] we show the 
same behavior for the plane wave formulation as in Fig. 
[2] for the density formulation, where we demonstrate that 
the total free energy is a constant of motion. Figure [8] 
shows the error as measured by the amplitude of the total 
energy vs. the free energy fluctuations for various lengths 
of the integration time step 5t. Only the free energy has 
the correct scaling with respect to the time step. 
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FIG. 5: (Color online) The fluctuations in the total free 
energy, E^p^ ^ Eyl + U — T^S, for a small Li cluster at Te = 
2, 000 K simulated using FreeON within the density matrix 
extension. The Pulay force has been calculated either with 
the original zero-temperature expression, Eq. (47 1, or its finite 
temperature generalization, Eq. (29 1. 



C. Fermi operator expansion 

The temperature-dependent density matrices for the 
simulations in Fig. |2] were constructed using the "exact" 
analytic Fermi-Dirac distribution of the states. If we in- 
stead use the recursive Fermi operator expansion of Al- 
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FIG. 6: (Color online) The change in total free energy as 
a function of time for a periodic 8 atom silicon cell for con- 
ventional Born-Oppenheimer MD as implemented in VASP 
and extended Lagrangian free energy MD. Both simulations 
using a time step of 0.25 fs, and a energy convergence cri- 
teria of 0.5 meV. The conventional Born-Oppenheimer MD 
show a clear unphysical systematic drift in comparison to the 
extended Lagrangian free energy MD. 
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FIG. 7: (Color online) The fluctuations for the total nuclear 
and potential energy, Ek + U, in comparison to the total free 
energy, Ek + U + TeS, of a periodic 8 atom silicon system. 
The simulation was carried out using extend Lagrangian free 
energy MD in VASP. Using a time step of 0.25 fs and en- 
ergy convergence criteria of 5 x 10~^° eV with ultrasoft pseu- 
dopotentials. Both the nuclear and electronic had an average 
temperature of about 500 K. 



gorithm [T] we find no significant deviation compared to 
the Fermi-Dirac result when we use more than 5 steps in 
the recursion. Figure |9] shows the free energy calculated 
with the exact result vs. the recursive expansion. The 
example corresponds to the simulation of the lower panel 
in Fig. [2] Using 5 recursion steps gives only a small con- 
stant shift in the free energy. The accuracy for a specific 
expansion order m will depend on the temperature Tg. In 
Fig. [To] we show the error in the entropy as a function of 
m for three different temperatures for a small Li cluster. 



FIG. 8: (Color online) The error amplitudes as a function of 
time step for the total nuclear and potential energy, Ek + U, 
and the total free energy, Ek + U + TeS, of a 8 atom periodic 
silicon cell 
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FIG. 9: (Color online) The fluctuations in the total free 
energy, _E|?' = Ek + U ~ T^S, for a water molecule, H2O 
(PBE0/3-21G), using either the exact Fermi-Dirac distribu- 
tion or the recursive Fermi operator expansion. Algorithm [l] 
with 5 (m = 5) or 6 (m — 6) recursion steps. The simula- 
tion is based on a hybrid density functional (PBEO [5H]) for 
Te = 10, 000 K and Tion ~ 300 K using 2 SCF cycles per time 
step, 5t — 0.5 fs, within the density matrix extension. 



In this case we find that the relative accuracy is reduced 
for lower temperatures, though the rate of convergence 
as a function of the number of expansion steps m is the 
same. 



VI. SUMMARY AND DISCUSSION 

Free energy generalizations of regular Born- 
Oppenheimer molecular dynamics for simulations 
at finite electronic temperatures are currently widely 
used [m [121 [M] . In this paper first principles free 



10 



Appendix A: Recursive Fermi operator expansion 
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FIG. 10: (Color online) The relative error in the entropy as a 
function of the number of expansion steps m in the recursive 
Fermi operator expansion, Algorithm [l] for three different 
electronic temperatures Tg. 



There are at least 19 different (dubious) ways to cal- 
culate matrix exponentials |65j . Consider 



which after a first order Taylor expansion gives 
— lim 



ca \ 2n — X 



(Al) 



(A2) 



This means that the Fermi-Dirac distribution function, 
$(a;), can be approximated as 

•^{x) = [e^ + l]-i = lim ^^"^ ~ ^j" (A3) 

' ^ J «^oo (2ri + x)" + (2n - a;)" ^ ' 

such that 



energy molecular dynamics schemes were developed 
based on extended Lagrangian Born-Oppenheimer 
molecular dynamics [21 SI UHl IS] . The formulation was 
given both for density matrix and plane wave represen- 
tations. The density matrix formulation has numerically 
convenient expressions for the electronic forces as well as 
a recursive Fermi operator expansion that are well suited 
for reduced complexity calculations. Our extended 
Lagrangian free energy molecular dynamics enables 
efficient and accurate molecular dynamics simulations 
based on Hartree-Fock and density functional theory (or 
any of their many extensions) both for localized orbital 
and plane wave pseudopotential calculations with a 
temperature-dependent electronic degrees of freedom 
and nuclear trajectories that evolve on the self-consistent 
free energy potential surface. 
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<^{2n - -inx) 



a;" -f (1 - xy 
for large values of n. The function 

fn{x) 



(1 - X)" 

has the important recursive property. 



(A4) 



(A5) 



(A6) 



which enables a rapid high-order expansions. The follow- 
ing recursive approximation can therefore be used: 

$[/3(£, -m)] = *(2n-4na;,) 

«/«(x.) = /2(/2(.../2(a:,)...)), (A7) 



where 



(A8) 



and with the recursion repeated m times, i. e. for n = 2™ 
The density matrix at finite electronic temperatures, 



D = 



= $ [P{H - fil)] , (A9) 



can thus be calculated using the recursive grand canoni- 
cal Fermi operator expansion, 

= /2 . . /2 (/2 Q/ - 2-(2+™)/3(iJ - A^/))) . . , 

(AlO) 

which is generated by Algorithm [l] 
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Appendix B: Constant energy integration: a 
practical guide 

1. Underlying approximations 



frequency and low- frequency modes |66| . The problem 
is particularly significant when there is a large difference 
between the atomic masses [57] , 

3. Practical guide 



Even if a correct geometric and long-term energy con- 
serving integration scheme, such as the Verlet algorithm, 
is used both for the nuclear and the electronic degrees 
of freedom, a systematic energy drift may still appear in 
a simulation. The most probable cause is the approxi- 
mate formulation of the underlying dynamics. Any self- 
consistent first principles calculation involves many dif- 
ferent numerical approximations, of which some can be 
physically motivated, others not. In the ideal case we 
would perform all calculations numerically exact, but for 
an approximate underlying dynamics that is sufficiently 
close to our effective single particle theory. The dominat- 
ing major approximations due to incomplete SCF conver- 
gence and the finite time steps are treated within the ex- 
tended Lagrangian framework and the geometric integra- 
tion presented here. Other approximations, such as trun- 
cated multipole expansions, finite arithmetics, threshold- 
ing of small matrix elements, and finite grid errors, have 
not been considered. First principles electronic structure 
codes have in general never been able to generate exact 
energy conserving molecular dynamics. To obtain a well 
working scheme that is stable and energy conserving un- 
der approximate SCF convergence may thus involve ex- 
tensive prior testing, tuning and debugging, even if the 
dissipative electronic force discussed above is applied. In 
this way, an energy conserving extended Lagrangian free 
energy (or Born-Oppenheimer) molecular dynamics sim- 
ulation is an excellent test and signature of a well working 
code. 



2. Thermostats 

In regular Born-Oppenheimer molecular dynamics sim- 
ulations, various forms of thermostats are applied to gen- 
erate, for example, stable constant temperature (NVT) 
or constant pressure (NVP) ensembles. In this case the 
problem with a systematic energy drift may not be seen 
directly, but the underlying dynamics is, of course, still 
wrong, even if a hypothetical "exact" thermostat could 
be used. By using an ever increasing number of SCF cy- 
cles the underlying energy drift can be kept low, but this 
is computationally expensive. Thermostats typically act 
through a rescaling of the forces or the velocities of the 
nuclear degrees of freedom. Unfortunately, this approach 
may lead to an unphysical energy transfer between high- 



Depending on the accuracy of the electronic structure 
code used in the molecular dynamics simulation, we give 
the following practical guide for a constant free energy 
simulation: 1) Use the extended Lagrangian formulation 
[3] and a geometric integration scheme [H [25H?7| for the 
integration of both the nuclear and the electronic degrees 
of freedom. 2) If the calculation is noisy, use electronic 
dissipation [34], e.g. as in Eq. (35 1. 3) If the simula- 



tion still shows excessive fluctuations in the energy or a 
systematic drift we suggest a careful review and tuning 
of the numerical approximations in the underlying elec- 
tronic description. 4) If the energetics still is unstable, 
and then only as a last resort, one could use a stabiliz- 
ing velocity or force rescaling. In that case we suggest 
the following velocity rescaling, which is constructed to 
minimize the perturbation of the molecular trajectories: 



EK{n) 



EkU) 



(Bl) 

Vi -> c„Vi. (B2) 

Here is the velocity of particle i, which is rescaled 
by the constant c„ at time step n. Ep is the desired 
target (free) energy, E^p*-{n) — Ek + U — TgS is the total 
(free) energy and EK{n) is the nuclear kinetic energy at 
time step n. The constant r is a chosen relaxation time 
and 7 is a weight factor for the integrated error term. 
The rescaling above is loosely based on the Berendsen 
thermostat f68 including an integrated error term (if we 
chose I7I > 0) as is often used in conventional control 
theory. The velocity rescaling factor c„ is typically small 
and if the Verlet integration scheme is used, c„ scales with 
the length of the integration step St as c„ = l + 0{5t^) or 
Cn = 1 + 0{St'^), depending on how the relaxation time 
T is chosen. For higher-order integration schemes, c„ = 
l+0{St'') for larger integers of k. This means that even if 
the velocity rescaling is unphysical, the error can at least 
be kept very small, with only minor perturbations of the 
molecular trajectories. The velocity rescaling above was 
not necessary for any of the simulations performed for 
this article. 
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